
module mglobal_smith

        use mConstants
        use mObject
        use mUniverse
        use mLPS
        use mPSFittingFunctions
        use mMisc

        type(tUniverse) :: xUniverse
        type(tLPS) :: xLPS
        type(tSmith) :: xSmith

end module

subroutine def_universe(om_m,om_v,om_b,h0,As,n,w,wa)

        use mConstants
        use mObject
        use mUniverse
        use mLPS
        use mPSFittingFunctions
        use mMisc
        use mglobal_smith

        implicit none

        real(D), intent(in) :: om_m,om_v,om_b,h0,As,n,w,wa

        !print *, 'start def_universe:'
        !print *, 'om_m = ',om_m
        !print *, 'om_v = ',om_v
        !print *, 'om_b = ',om_b
        !print *, 'h0 = ',h0
        !print *, 'As = ',As
        !print *, 'n = ',n
        !print *, 'w = ',w
        !print *, 'wa = ',wa

        call new(xUniverse, omega_m0=om_m, omega_v0=om_v, omega_b0=om_b, &
                h=h0, A_s=As, w_0=w, w_a=wa)
        !print *, 'after new universe'
        call new(xLPS, universe = xUniverse, spectralIndex=n, normalization = 'CB') 
        !print *, 'after new lps'
        call new(xSmith, xLPS)
        !print *, 'after new smith'

end subroutine def_universe

real(D) function smithp3d(k,z)

        use mConstants
        use mObject
        use mUniverse
        use mLPS
        use mPSFittingFunctions
        use mMisc
        use mglobal_smith

        implicit none

        real(D), intent(in) :: k,z

        real(D) :: a

        a = 1.0_D/(1.0_D+z)

        !print *, 'start smithp3d: k = ',k,' a = ',a

        call smithEtAl( xSmith, k, a, smithp3d )

        !print *, 'done smithp3d: smithp3d = ',smithp3d

end function smithp3d

real(D) function linearp3d(k,z)

        use mConstants
        use mObject
        use mUniverse
        use mLPS
        use mPSFittingFunctions
        use mMisc
        use mglobal_smith

        implicit none

        real(D), intent(in) :: k,z

        real(D) :: a, growth

        a = 1.0_D/(1.0_D+z)

        !print *, 'start linearp3d: k = ',k,' a = ',a

!        linearp3d = P_L( xLPS, k, a )
        growth = growthFactor( xUniverse, a )
        linearp3d = growth * growth * P_L( xLPS, k, a)

        !print *, 'done linearp3d: linearp3d = ',linearp3d

end function linearp3d

real(D) function pdp3d(k,z)

        use mConstants
        use mObject
        use mUniverse
        use mLPS
        use mPSFittingFunctions
        use mMisc
        use mglobal_smith

        implicit none

        real(D), intent(in) :: k,z

        real(D) :: a

        a = 1.0_D/(1.0_D+z)

        !print *, 'start pd: k = ',k,' a = ',a

        pdp3d = peacockDodds( xLPS, k, a )

        !print *, 'done pd: pd = ',pdp3d
       
end function pdp3d

